Bounds in Probability

Adapted from: SOSTOOLS' SOSDEMO8 (See Section 4.8 of SOSTOOLS User's Manual)

The probability adds up to one.

μ0 = 1

The mean is one.

μ1  = 1

The standard deviation is 1/2.

σ = 1/2

0.5

The second moment E(x^2) is:

μ2 = σ^2 + μ1^2

1.25

We define the moments as follows:

using DynamicPolynomials
@polyvar x
monos = [1, x, x^2]
using SumOfSquares
μ = measure([μ0, μ1, μ2], monos)

Measure{Float64, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}([1.0, 1.0, 1.25], DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}[1, x, x²])

We need to pick an SDP solver, see here for a list of the available choices. We use SOSModel instead of Model to be able to use the >= syntax for Sum-of-Squares constraints.

using CSDP
solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)
model = SOSModel(solver);

We create a polynomial with the monomials in monos and JuMP decision variables as coefficients as follows:

@variable(model, poly, Poly(monos))

\[ ({\_}_{1}) + ({\_}_{2})x + ({\_}_{3})x^{2} \]

Nonnegative on the support:

K = @set 0 <= x && x <= 5
con_ref = @constraint(model, poly >= 0, domain = K)

\[ ({\_}_{1}) + ({\_}_{2})x + ({\_}_{3})x^{2} \text{ is SOS} \]

Greater than one on the event:

@constraint(model, poly >= 1, domain = (@set 4 <= x && x <= 5))

\[ ({\_}_{1} - 1) + ({\_}_{2})x + ({\_}_{3})x^{2} \text{ is SOS} \]

The bound (we use LinearAlgebra for the ⋅ syntax for the scalar product):

using LinearAlgebra
@objective(model, Min, poly ⋅ μ)

$ {_}{1} + {\}{2} + 1.25 {\}_{3} $

We verify that we found a feasible solution:

optimize!(model)
primal_status(model)

FEASIBLE_POINT::ResultStatusCode = 1

The objective value is 1/37:

objective_value(model)

0.027027027945194515

The solution is (12x-11)^2 / 37^2:

value(poly) * 37^2

\[ 120.9999573610103 - 263.99994311062636x + 143.99998960526992x^{2} \]

This page was generated using Literate.jl.